// retrieval.h

#ifndef DODO_FERNALD_H
#define DODO_FERNALD_H

#include <cmath>
#include <vector>
#include <stdexcept>

#include <quigon.h>

namespace quigon
{

// required fields in para:
//	Lidar_Profile			"distance"
//	Lidar_Cross_Ref			"distance_i"
//	Link_Stuff<Time_Serie>	"lidar_constant", "lidar_ratio", "sin_elev_angle"
//	size_t					"zero_meter_index", "first_usable_index"
//	double					"bin_size"
template<class Tline, class Tpara>
class Fernald
{
public:
	Fernald(
	        const Tpara& para0,
	        const Atmos_Profile& beta_m0,
	        double max_distance0 = 10.0,
	        double alt_sa0 = -999.9 , // minus mean not using it.
	        double alt_lidar_constant0 = -999.9,
	        bool kick_minus0 = false
	) :
			beta_m(beta_m0),
			max_distance(max_distance0),
			alt_sa(alt_sa0),
			alt_lidar_constant(alt_lidar_constant0),
			current_sa(-999.9),
			exp_term(para0.size),
			old_x_exp(para0.size),
			cal(para0.size),
			kick_minus(kick_minus0),
			para(para0)

	{
		acc_beta_m = gen_acc_profile(beta_m, para.zero_meter_index, para.bin_size);
	}

	double& get_sa() { return alt_sa; }
	double& get_lc() { return alt_lidar_constant; }
	const double& get_sa() const { return alt_sa; }
	const double& get_lc() const { return alt_lidar_constant; }
	bool operator()(Tline& line)
	{
		try
		{
			if (!line.is_valid())
				return true;

			if ((line.status() & Status::EX))     // when it's already extinction line;
			{
				return true;
			}

			size_t max_index = para.distance_i[max_distance];
#ifdef VERBOSE
			std::cout << "max_index : " << max_index << std::endl;
#endif
			if (max_index >= line.size())
				max_index = line.size() - 1;

			// caution about max_index == 0;
			// should add a judge here
			line.resize(max_index + 1);

			std::vector<typename Tline::value_type> old(line.begin(), line.end());    // have to have a temp data. so copy a minimum one.

			double sa = para.lidar_ratio()[line.t()];  // lidar_ratio is a Link_Stuff<Time_Serie>
			double lc = para.lidar_constant()[line.t()];  //Link_Stuff<Time_Serie>

			// determine whether use local sa and lc.
			if (alt_sa > 0.0)
				sa = alt_sa;
			if (alt_lidar_constant > 0.0)
				lc = alt_lidar_constant;

			double ce_term = 1.0;
			if (!(line.status() & Status::LC))
				ce_term *= lc;
			if (!(line.status() & Status::ENERGY))
				ce_term *= line["energy"];

			size_t start_index = para.zero_meter_index;
			double bin_size = para.bin_size; // make the following more clear

			if (current_sa != sa)    // refresh exp_term;
			{
				current_sa = sa;
				for (size_t i = start_index; i <= max_index; ++i)
					exp_term[i] = exp(-2.0 * (sa - SM) * acc_beta_m[i]);
			}

			for (size_t i = start_index; i <= max_index; ++i)
			{
				old_x_exp[i] = old[i] * exp_term[i];
				if (i == start_index)
					cal[i] = old_x_exp[i] * bin_size * 0.5;
				else
					cal[i] = cal[i - 1] + (old_x_exp[i - 1] + old_x_exp[i]) * bin_size * 0.5;

//                    for ( size_t j = start_index; j <= i; ++j )
//                    {
//                        cal += old[j] * exp_term[j];
//                    }
//                    cal *= bin_size;

				line[i] = (old_x_exp[i]
				           / (ce_term - 2.0 * sa * cal[i]) - beta_m[i])
				          * sa;

				if (kick_minus && (line[i] < 0.0))
					line[i] = 0.0;
			}

			line.status() = Status::EX;
			return true;
		}
		catch (...)
		{
			func_err("Fernald::operator()");
			throw;
		}
	}

private:
	double max_distance;
	double alt_sa;  // here we'll call lidar_ratio sa for short;
	double alt_lidar_constant;

	const Atmos_Profile& beta_m;
	Atmos_Profile acc_beta_m;

	double current_sa;
	std::vector<double> exp_term;
	// exp( -2.0 * ( sa - SM ) * acc_beta_m[i] ); To make it faster !
	std::vector<double> old_x_exp;
	std::vector<double> cal;

	bool kick_minus;

	// from para;
	const Tpara& para;
};

template<class Tline, class Tpara>
class Fernald_Priori
{
public:
	Fernald_Priori(
	        const Tpara& para0,
	        const Atmos_Profile& beta_m0,
	        const Time_Serie& priori_distance0,
	        const Time_Serie& priori_value0,
	        double max_distance0 = 10.0,
	        double alt_sa0 = -999.9 , // minus mean not using it.
	        double alt_lidar_constant0 = -999.9,
	        bool kick_minus0 = false
	) :
			beta_m(beta_m0),
			max_distance(max_distance0),
			alt_sa(alt_sa0),
			alt_lidar_constant(alt_lidar_constant0),
			current_sa(-999.9),
			exp_term(para0.size),
			old_x_exp(para0.size),
			cal(para0.size),
			kick_minus(kick_minus0),
			para(para0),
			priori_distance(priori_distance0),
			priori_value(priori_value0),
			current_priori_index(-999)
	{
	}

	double& get_sa() { return alt_sa; }
	double& get_lc() { return alt_lidar_constant; }
	const double& get_sa() const { return alt_sa; }
	const double& get_lc() const { return alt_lidar_constant; }

	bool operator()(Tline& line)
	{
		try
		{
			if (!line.is_valid())
				return true;

			if ((line.status() & Status::EX))     // when it's already extinction line;
			{
				return true;
			}

			size_t max_index = para.distance_i[max_distance];
			size_t max_index_mem = max_index;
			size_t priori_index = para.distance_i[priori_distance[line.t()]];
#ifdef VERBOSE
			std::cout << "max_index : " << max_index << std::endl;
			std::cout << "priori_index : " << priori_index << std::endl;
#endif
			if (priori_index > line.size() - 1)
			{
				func_err("Fernald_Priori::()", "priori is too far");
				return false;
			}

			if (max_index < priori_index)
				max_index = priori_index;

			// caution about priori_index == 0;
			// should add a judge here
			line.resize(max_index + 1);

			std::vector<typename Tline::value_type> old(line.begin(), line.end());    // have to have a temp data. so copy a minimum one.

			double sa = para.lidar_ratio()[line.t()];  // lidar_ratio is a Link_Stuff<Time_Serie>
			double lc = para.lidar_constant()[line.t()];  //Link_Stuff<Time_Serie>

			// determine whether use local sa and lc.
			if (alt_sa > 0.0)
				sa = alt_sa;
			if (alt_lidar_constant > 0.0)
				lc = alt_lidar_constant;

			size_t start_index = para.zero_meter_index;
			if (start_index > priori_index)
			{
				func_err("Fernald_Priori::()", "priori under ground?");
				return false;
			}

			double bin_size = para.bin_size; // make the following more clear
			double priori_value_now = priori_value[line.t()];

			if (current_priori_index != priori_index)   // refresh acc_beta_m.
			{
				current_priori_index = priori_index;
				acc_beta_m = gen_acc_priori_profile(beta_m, priori_index, start_index, line.size() - 1, para.bin_size);
			}
			if (current_sa != sa)    // refresh exp_term;
			{
				current_sa = sa;
				for (size_t i = start_index; i < line.size(); ++i)
					exp_term[i] = exp(-2.0 * (sa - SM) * acc_beta_m[i]);
			}
			for (size_t i = priori_index; i >= start_index; --i)
			{
				old_x_exp[i] = old[i] * exp_term[i];
				if (i == priori_index)
					cal[i] = 0.0;
				else
					cal[i] = cal[i + 1] - (old_x_exp[i + 1] + old_x_exp[i]) * bin_size * 0.5;

//                    for ( size_t j = start_index; j <= i; ++j )
//                    {
//                        cal += old[j] * exp_term[j];
//                    }
//                    cal *= bin_size;

				line[i] = (old_x_exp[i]
				           /	(
				                   old[priori_index] / (beta_m[priori_index] + priori_value_now / sa)
				                   - 2.0 * sa * cal[i]
				           )
				           - beta_m[i])
				          * sa;

				if (kick_minus && (line[i] < 0.0))
					line[i] = 0.0;
			}
			for (size_t i = priori_index + 1; i < line.size(); ++i)
			{
				old_x_exp[i] = old[i] * exp_term[i];
				if (i == priori_index)
					cal[i] = 0.0;
				else
					cal[i] = cal[i - 1] + (old_x_exp[i - 1] + old_x_exp[i]) * bin_size * 0.5;

//                    for ( size_t j = start_index; j <= i; ++j )
//                    {
//                        cal += old[j] * exp_term[j];
//                    }
//                    cal *= bin_size;

				line[i] = (old_x_exp[i]
				           /	(
				                   old[priori_index] / (beta_m[priori_index] + priori_value_now / sa)
				                   - 2.0 * sa * cal[i]
				           )
				           - beta_m[i])
				          * sa;

				if (kick_minus && (line[i] < 0.0))
					line[i] = 0.0;
			}

			line.resize(max_index_mem + 1);
			line.status() = Status::EX;
			return true;
		}
		catch (...)
		{
			func_err("Fernald::operator()");
			throw;
		}
	}

private:
	double alt_sa;
	double alt_lidar_constant;
	double max_distance;
	const Time_Serie& priori_distance;
	const Time_Serie& priori_value;
	const Atmos_Profile& beta_m;
	Atmos_Profile acc_beta_m;

	double current_sa;
	size_t current_priori_index;
	std::vector<double> exp_term;
	std::vector<double> old_x_exp;
	std::vector<double> cal;

	bool kick_minus;

	const Tpara& para;
};

}


#endif

